#!/usr/bin/python
from zplot import *
import xml.etree.ElementTree as ET
import parse
from math import *
import pylab as pl
import random

def geo_scale(lat, lon):
	M = 6367450.
	phi = lon/360.*pi*2
	psi = lat/360.*pi*2
	lat_meters = 111133-560*cos(2*phi)+1.2*cos(4*phi)
	lon_meters = pi*M*cos(psi)/180.
	ratio = lat_meters/lon_meters
	return (lat_meters, lon_meters)

#call geo_scale first and then let
#let x = lat_meters*lat and y = lon_meters*lon

def total_curvature(x, y):

    n = len(x)

    total_k = 0
    length = (x[0]**2+y[0]**2)**0.5+(x[n-1]**2+y[n-1]**2)**0.5

    for i in range(1, n-1):

        t1x = x[i+1]-x[i]
        t1y = y[i+1]-y[i]

        t0x = x[i]-x[i-1]
        t0y = y[i]-y[i-1]

        s1 = (t1x**2+t1y**2)**0.5
        s0 = (t0x**2+t0y**2)**0.5

        if s1 == 0 or s0 == 0:
            continue

        T1x = t1x/s1
        T1y = t1y/s1

        T0x = t0x/s0
        T0y = t0y/s0

        s = 0.5 * ((x[i+1]-x[i-1])**2+(y[i+1]-y[i-1])**2)**0.5

        kx = (T1x-T0x)/s
        ky = (T1y-T0y)/s

        k = (kx**2+ky**2)**0.5*s

        #k = (kx**2+ky**2)*s

        total_k += k
	return total_k

def length(x,y):
    n = len(x)
    L = 0
    for i in range(1, n):
        L += ((x[i]-x[i-1])**2+(y[i]-y[i-1])**2)**0.5
    return L


class MapVis(object):
    def __init__(self, path):
        self.data = parse.Parse(path)

    def average_curvature(self):

        L = 0
        curvature = 0

        curvatures = []
        street_lengths = []

        n = len(self.data.ways)
        w1 = self.data.ways[0]
        lat_meters, lon_meters = geo_scale(w1.nodes[0][0], w1.nodes[0][1])
        #lat_meters, lon_meters = 1,1
        print n
        sumx = 0
        sumy = 0
        for k in range(0, min(n, n)):
            w = self.data.ways[k]
            x = []
            y = []
            if len(w.nodes) > 3:
                startx=w.nodes[0][0]
                starty=w.nodes[0][1]

                m = len(w.nodes)

                endx=w.nodes[m-1][0]
                endy=w.nodes[m-1][1]

  
                for i in range(len(w.nodes)):  
                    n1 = w.nodes[i]
                    #x.append(n1[0]*lat_meters)
                    #y.append(n1[1]*lon_meters)
                    x.append((n1[0]-startx)*lat_meters)
                    y.append((n1[1]-starty)*lon_meters)
                    pl.plot(x,y)

                pl.plot(x,y)

                sumx += (w.nodes[len(w.nodes)-1][0] - startx)
                sumy += (w.nodes[len(w.nodes)-1][1] - starty)

                k = total_curvature(x,y)
                curvatures.append(log(k/length(x,y)+10**(-10))/log(10))
                #curvatures.append(k/length(x,y))
                curvature += k
                street_length=length(x,y)
                L += street_length
                street_lengths.append(street_length)
                #print street_length

       
        print "average curvature: \t", curvature/L
        print "geometric mean of averages:\t", 10**(sum(curvatures)/len(curvatures))
        pl.hist(curvatures, 20)
        print len(curvatures)

        pl.show()
        



def main():
    filename = "data/wichita.xml"
    print filename
    mv = MapVis(filename)
    mv.average_curvature()


if __name__ == '__main__':
    main()
